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ABSTRACT 

Gravitational instabilities play an important role in galaxy evolution and in shaping the inter¬ 
stellar medium (ISM). The ISM is observed to be highly turbulent, meaning that observables 
like the gas surface density and velocity dispersion depend on the size of the region over 
which they are measured. In this work we investigate, using simulations of Milky Way-like 
disc galaxies with a resolution of ~ 9 pc, the nature of turbulence in the ISM and how this af¬ 
fects the gravitational stability of galaxies. By accounting for the measured average turbulent 
scalings of the density and velocity fields in the stability analysis, we can more robustly char¬ 
acterize the average level of stability of the galaxies as a function of scale, and in a straightfor¬ 
ward manner identify scales prone to fragmentation. Furthermore, we find that the stability of 
a disc with feedback-driven turbulence can be well described by a “Toomre-like” Q stability 
criterion on all scales, whereas the classical Q can formally lose its meaning on small scales 
if violent disc instabilities occur in models lacking pressure support from stellar feedback. 

Key words: instabilities - turbulence - ISM; general - ISM: kinematics and dynamics - 
ISM: structure - galaxies: ISM 


1 INTRODUCTION 


Today, over three decades after the pioneering work by |Larson| 
( |1981^ , observations and simulations of the interstellar medium 
(ISM) are revealing its turbulent nature with higher and higher fi¬ 
delity (see review by e.g. |Mac Low & Klessen||200 4{ [Elmegre^ 
|& Scalo|2004| ). Turbulence is not only thought to play an impor¬ 
tant role in controlling star formation in molecular clouds jMcKee 


& Ostriker|2007| [Padoan & Nordlund|2011||Federrath & Klessm 


2012| l, but also on galactic scales (|Renaud^^e^^alJj201A|and in 
shaping the interstellar medium (ISf^ jStanimi rovic et ^|1999[ 

lElmegreen et al.|2001[[Boumaud et al.|2010| l. The importance of 

galactic scale turbulence has in the past decade also been revealed 
in the early Universe; in gas rich high redshift galaxies, the ob¬ 
served levels of gas turbulence are much higher than in the lo¬ 
cal Universe (Shapiro et al.||200^ |Edrster Schreiber et al.||200^ 


ISwinbank et ah |2011^ , explaining the observed ubiquity of super 

star forming clumps (Elmegreen et al.|2009[ [Agertz et al. 


2009b |Bournaud & Elmegreen||2009[ Dekel et al. |2009 Genzel 


etal.| 26 ll|l. 

One of many fundamental aspects of ISM turbulence is the 
existence of scaling relations between observables, such as the col¬ 
umn density (E), the ID velocity dispersion (cr), and the size of the 
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region (t) over which such quantities are measured: 

E oc (T oc t. (1) 

The values of exponents a and b depend on which ISM component 
and the range of scales that are considered. In this work we focus 
mainly on the cold neutral gas; neutral hydrogen (HI), and molec¬ 
ular gas, dominated by molecular hydrogen (H 2 , observed via the 
tracer molecule CO), which is known to be supersonically turbu¬ 
lent and plays an important role in the gravitational instability of 
galactic discs (e.g.|Lin & Shu|1966[ |Jog & Solomon|1984[|Bertin| 
|& Romeo|1988| >. 

In molecular gas, the scaling exponents are a « 0 and t> ~ |, 
up to scales of several 100 p<0 This pair of exponents are often 
referred to as “Larson’s scaling laws” after the discovery by |Lar- | 
|sonH198T) (see also [Solomon et aT||1987| >. In fact, both Galactic 
and extragalactic Giant Molecular Clouds (GMCs) are fairly well 
described by Larson’s scaling laws, although with large uncertain- 


ties (e.g. Bolatto et al. 2008 Heyer et al.|2009 

|Hughes et al. 

2010 

Kauffmann et al.|2010l|Lombardi et al.|2010l 

Sanchez et al. 

2010 

Roman-Duval et al.||2010l|Ballesteros-Paredes et al.||2011l 

Beau- 

mont et al.|2012|>, as well as in the dense star-forming clumps in 


high redshift galaxies jSwinbank et al.|201 Ijl. 


Simulations of GMCs forming in galactic discs have recently 


^ Note that a = 0 is expected for isolated clouds in gravitational equi¬ 
librium, as the cloud mass M oc together with a oc ® gives 
Yi Mj= constant. 
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started producing scaling relations compatible with Larson’s re¬ 
lations (e.g. [Hopkins et al.]|2012[ |Dobbs|2014), but see |Fujimoto| 


et al. j2014^ for models predicting steeper relations (E(£), (t(£) oc 


£). State-of-the-art numeric al simulations of supersonic turbulence 
suggest a ~ and fo ~ ^ ( Fleck] 1996 ||Kowal & Lazarianl2007| 


Kritsuk et al.|2007l ISchmidt et al.|2008||Price & Federrath|20l^ 


Kritsuk et al.|2013] l. Other recent numerical work suggest that the 
scaling exponent a may be significantly affected by the nature of 
the turbulence forcing ( Federrath et al.|20Tol >, magnetic fields, and 
self-gravity ^Collins et al. 2012] , 

In Hi, the scaling exponents are a ~ | and 6 ~ | up to 
several kpc ( |Roy et al.|2008|>, altho ugh, as for molecular gas, with 
large uncertainties ( Kim et al.|2007| >. Observed power spectra of HI 
intensity fluctuations in nearby galaxies are compatible with a Kol¬ 
mogorov scaling for both a and E (e.g. Stanimirovice£^J1^9^ 


Lazarian & Pogosyan|2000{ [Elmegreen et al.|2001[ Begum et al.| 


2006[|Bournaud et al. [20lo| [Zhang et al.|2012[|Dutta et al.|2013| >, 

with similar results for other ISM components on large scales (e.g. 
dust and CO in M33, [Comhes et al.[2012) . Furthermore, HI power 
spectra are often found to be shallower on large scales, with a 
break around i ~ the disc scale height, possibly indicating a transi¬ 
tion from 3D turbulence on small scales to 2D turbulence on large 
( [Dutta et al.|2008||2009[|Block et al.|2010[[Boumaud et al.|2010] . 

While the turbulent nature of the ISM is well established, it 
is rarely accounted for in theoretical works when evaluating the 
gravitational stability of galactic discs. Instead, cr and E are asso¬ 
ciated with smoothed quantities on galactic scales (~ kilo-parsecs) 
(hut see [Elmegreen|[19 96[ [Begelman & Shlos man[[200^ [H opkins[ 
[2012[ >. A few analytical studies investigating the effect of ISM tur¬ 
bulence on gravitational stability have been carried out. Romeo et 
al. (2010) explored a range of values for a and h and showed that 
turbulence has an important effect on the gravitational instability 
of the disc; it excites a rich variety of stability regimes, several of 
which have no classical counterpart. Followup studies b 3 HHoffinam 


& Romeo| ( [2012] and [Romeo & Falstad[j2013| > (see also [Shadmehri 


& Khajenabi 2012[> extended this framework to turbulent multi- 


component (gas-l-stars) discs and applied it to the THINGS galaxy 
sample ^Walter et al.[2008T l. Their analysis showed that H 2 plays a 
significant role in disc (in)stability by dynamically decoupling and 
dominating the onset for gravitational instability even at distances 
as large as half the optical radius. 

The goal of this work is to in greater details explore the inter¬ 
play between gas turbulence and disc stability by extending previ¬ 
ous work in a number of ways: 


(i) We perform numerical simulations of a Milky Way-like 
galaxy with a resolution of ~ 9 pc, where we model the same 
galaxy in two markedly different ways: (1) without any stellar feed¬ 
back, leading to rapid gas fragmentation into a population of star 
forming giant molecular clouds (GMCs), and (2) with efficient stel¬ 
lar feedback which acts to disperse the GMCs, drives interstellar 
turbulence, and regulates the rate of star formation. These two mod¬ 
els show two extremes of galaxy evolution, and are useful platforms 
on which to understand, and characterize, the role of turbulence in 
realistic disc galaxies. 

(ii) We apply a classical stability analysis on kpc-scales to both 
models and demonstrate how it fails to qualitatively separate the 
two systems, despite the dramatically different ISM morphology. 
We show that the this is in agreement with observations of nearby 
spiral galaxies. 

(iii) By accounting for the scale-dependent nature of E and a in 
a multi-component stability analysis, we demonstrate how we can 


more robustly characterize the disc stability and dynamical state of 
the galaxies. 

The rest of the paper is organized as follows. The classical 
framework for understanding the stability of two-component thick 
discs is outlined in Section l2.1l and lT2l In Section|23]we describe 
the numerical hydro-t-A^-body method used for the galaxy simula¬ 
tions. In Section ED we perform a classical analysis followed by 
an accounting of turbulent scaling relations in Section [T2l In Sec¬ 
tion |3^ we discuss our results in a more general context of ob¬ 
served, and theoretically predicted, ISM scaling relations and what 
this means for understanding disc galaxies. Finally, we discuss and 
conclude our results in Section|4] 


2 METHOD 


2.1 Stability diagnostics 

Consider a gas disc of scale height h, and perturb it with axisym- 
metric waves of frequency w and wavenumber k. The response of 
the disc is described by the dispersion relation 


2 2 

uj = n 


27rGE k 
1 + kh 


I 2,2 

+ a k . 


( 2 ) 


where k is the epicyclic frequency, E is the surface density at equi¬ 
librium, and a is the sound speed (|Romeo[[ 1 992[ [ 1 994|( (see also 
[Vandervoort| 1970) . The three terms on the right side of Eq. [^rep¬ 
resent the contributions of rotation, self-gravity and pressure. Eor 
kh 1, Eq. (1) reduces to the usual dispersion relation for an 
infinitesimally t hin gas disc (for a derivation, see e.g. [Binney 
Tremaine]|2008 '. Such a disc is unstable if and only if oj <0, 


which is equivalent to Q < 1, where Q is defined by 




ttGE 


(3) 


This quantity was first derived by [Toomi^ ( [ 1 964) for a thin disc of 
stars (where cr —>• cr^, i.e. the stellar radial velocity dispersion, and 
the denominator has tt replaced by 3.36), and we hence refer to it as 
Toomre’s Q. Erom now on we denote Eq.j^for stars and gas as Q* 
and Qg respectively. Eor kh ';$> 1, one recovers the case of Jeans 
instability with rotation, since E/h « 2p. In other words, scales 
comparable to h mark the transition from 2D to 3D stability. 

Much work has gone into characterizing gravitational instabil¬ 
ities in multi-component systems jBertin & Romeo|1988[|^meo[ 


[T^ [Elmegreen] 1995 [[ Jog| 1 996[[Rafikov|200T) rSuch sys 


terns are always more unstable than each component considered 
separately, and the interplay between the different components de¬ 
pends on the ratio between their velocity dispersions and their 
surface densities. [Romeo & Wiege^ POI l[ l introduced a simple 
and accurate approximation for the two-component Q parameter, 
which takes into account the stabilizing effect of disc thickness and 
predicts whether the local stability level is dominated by stars or 
gas. [Romeo & Ealstadj ( [2013) generalized this approximation to 
discs made of several stellar and/or gaseous components, and to the 
whole range of velocity dispersion anisotropies observed in galac¬ 
tic discs. In the two-component case, the Q stability parameter is 
given by 
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if T’*Q* > TgQg , 

Qthick 

1 
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w = - 

C 

1 
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+ cr| ■ 

( / 


it ; 

T ^ < 

1 + 0.6 f 

//) 

for 0 < UzlerR < 0.5 , 


0 . 8 + 0.7 

(s) 

for 0.5 < crz/(T_R < 1. 


The thin-disc limit, Qthin, can be recovered by setting 
= 0 and hence T = 1. As for the single component case, 
the criterion for instability is Qthick < 1 and Qthin < 1 - 


2.2 Accounting for turbulent scalings 

As discussed in Sect. 1, observations and theory have revealed that 
the interstellar medium (ISM) is highly turbulent, and many prop¬ 
erties of the ISM depend on the physical scale on which they are 
measured. As argued by [Romeo et al.|p010[ l, turbulence can be 
accounted for in the stability diagnostics by generalizing the dis¬ 
persion relation in Eq.[^ to account for scale-dependencies of the 
velocity dispersion cr and surface density E, i.e. 

u? = -2nGT.{k)k + o-'^{k)k'^, (7) 

where E(fc) and (j{k) are the mass column density and the ID ve¬ 
locity dispersion (accounting both for the thermal and turbulent 
component) measured over a region of size f = 27r/fc. Observa- 


tions and theoretical studies (see e.g. 

Larson||1981||Elmegreen & 

Scalo|2004l|McKee & Ostriker|2007l 

Kritsuk et al.|2007l|Romeo 

et al.| 2010 l indicate a power-law behaviour in these quantities, mo- 


tivating the following parametrization: 


If the disc has volume density p and scale height h, then E « 2pl 
for t'fh and E « 2ph for t>h. The range tfh corresponds 
to the case of 3D turbulence (relevant for GMCs and HI on small 
scales), whereas the range corresponds to the case of 2D tur¬ 
bulence (relevant for HI on large scales). The quantity £o = 27r/feo 
is the typical scale at which E and a are observed and quan¬ 
tities like the Toomre parameter Q is computed, so that Qo = 
kcto/ttGEo. For example, at an angular resolution of 6 " achieved 
for the "The HI Nearby Galaxy Survey” (THINGS) galaxy sample 
( [Walter et al.|2008l l, analysis is carried out on the spatial smoothing 
scale £o ~ 0.5 — 1 kpc ( [Leroy et al.[2008[ (. 

As we are interested in understanding the influence of turbu¬ 
lence on a two component disc (gas and stars), we need to analyze 
the joint dispersion relation which can then be expressed in the 
following form (e.g. [Jog & Solomon||1984[ [Hoffmann & Romeo| 

[20T2I 1: 

- Mf) (cu" - Mi) = (Pi - M?) (Pi - Mi) , (9) 

where the (turbulent) dispersion relation for each component i has 
the usual form 

Mi = id — 2'KGT,i{k) k -|- (Ji{k) d, (10) 


and 

'P'i=d+af{k)d, (11) 

with i =gas, sta|^ Note that Mi(fc) is the one-component disper¬ 
sion relation for potential-density waves, while Vi{k) describes 
sound waves modified by rotation and turbulence. Furthermore, the 
right-hand side of Fq.|^ measures the strength of the gravitational 
coupling between the two components, as Mf (fc) — Vi{k) is the 
self-gravity term of component i. Fq.|^is quadratic in d, with two 
real roots: 

uil = ^[Mj + Ml±VAj, (12) 

A = (M? - Mi)" -f 4 (Pi - M?) (P| - Mi) . (13) 

This means that the dispersion relation has two branches that do 
not cross, uj'^(k) 7 ^ uj1_{k), except possibly as fe —>■ 0 or fe —>■ 00 . 
We focus only on uj 1_ (k) in this work, as this branch is the only one 
which can represent gravitational instability (aj^ (k) is manifestly 
> 0). Note that Lot (k) is always smaller than each component in¬ 
dividually. 

In Section [ 3 ^ and [33] we use the above framework, together 
with numerical simulations, to demonstrate how turbulence affects 
the shape of the dispersion relation, and hence the condition for 
local gravitational instability (cj^ < 0 ). 


2.3 Numerical technique 


In order to carry out hydro-(-A-body simulations of galactic discs, 
we use the Adaptive Mesh Refinement (AMR) code RAMSES 
( [Teyssier[2002^ . The fluid dynamics of the baryons is calculated us¬ 
ing a second-order unsplit Godunov method, while the collisionless 
dynamics of stellar and dark matter particles is evolved using the 
particle-mesh technique ( [Hockney & Fastwood[198T| l, with gravi¬ 
tational accelerations computed from the gravitational potential on 
the mesh. The gravitational potential is calculated by solving the 
Poisson equation using the multi-grid method ^Guillet & Teyssierj 
[2011[ | for all refinement levels. The equation of state of the fluid is 
that of an ideal mono-atomic gas with an adiabatic index 7 = 5/3. 

The code achieves high resolution in high density regions 
using adaptive mesh refinement, where the refinement strategy is 
based on a quasi-Lagrangian approach in which the number of col¬ 
lisionless particles per cell is kept approximately constant. This 
allows the local force softening to closely match the local mean 
interparticle separation, which suppress discreteness effects (e.g., 
[Romeo et al.[200^ . An analogous refinement criterion is also used 
for the gas. 

The star formation, cooling physics and stellar feedback 


model adopted in our simulations is described in detail in Agertz 


[et al. I 2013^ (identical to the ”A11” model) and[Agertz & Kravtsov 


(2014 1 , and we refer the reader to those papers for details. Briefly, 


several processes contribute to the stellar feedback budget, as stars 
inject energy, momentum, mass and heavy elements over time via 
SNII and SNIa explosions, stellar winds and radiation pressure 
into the surrounding gas. Metals injected by supemovae and stellar 
winds are advected as a passive scalar and are incorporated self- 
consistently in the cooling and heating routine. 


^ The dispersion relation of an TV-component turbulent disc is 
— P?)/(a;^ ~ 'Pf) = J- t’s inferred from equaiton 

22 of Rafikov (2001). 
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One aspect differs from previous work; at the current numer¬ 
ical resolution (Aat ~ 9 pc), we are not guaranteed to always re¬ 
solve the cooling radiu^rcooi, i-e. the SN bubble radius for which 
radiative losses are expected to be important for each discrete SN 
event, leading to an underestimation of the impact of SNe feed¬ 
back. Instead of remedying this issue by solving two separate en¬ 
ergy equations, one for the thermal energy and one for the feed¬ 
back energy as in|Agertz & Kravtsov| j2014f, we adopt the model 
recently suggested by Kim & Ostriker 1 201^ (see also |Martizzi| 
|et al.|[20T4{ |Gatto et A.||2014 Simpson et al'|[20T4) . Here a SN 
explosion resolved by at least three grid cells (rcooi > 3Aa:) is 
initialized in the energy conserving phase by injecting the relevant 
energy (10®^ erg per SN) into the nearest grid cell. If this criterion 
is not fulfilled, the SN is initialized in its momentum conserving 
phase, i.e. the momentum generated during the energy conserving 
Sedov-Taylor phase is injected into to the 26 cells surrounding a 
star particle. 

[Blondin et al.| ( [T998t calculated the transition time from the 
energy conserving phase to the phase of shell formation, at which 
the cooling time equals the age of the remnant (fcooi = fsN), to 
be « 2.9 X 10"* yrs, where no is the ambient density 

and i? 5 i the thermal energy in units of 10®^ ergs. At this time, the 
momentum of the expanding shell is approximately 


psT ~ 2.6 X 10® Mq kms 


(14) 


For reasonable values of ambient densities, this is ~ 10 times 
greater than the initial ejecta momentum. |Kim & Ostriker| \2QIA ) 
and |Martizzi et al.H2014^ have shown, using detailed simulations of 
SNe explosions, that Eq.[T^holds even for more realistic, clumpy, 
environments. We hence use this relation for the injected momen¬ 
tum per individual SN explosion when the cooling radius is not 
resolved. 


2.3.1 Simulations 


We model the non-linear evolution of an entire Milky Way-like 
galactic disc. This setup is now the standardized test for the AGORA 
code comparison project ^Kim et al.|20T4), and is a higher resolu¬ 
tion version of the galaxy analyzed in|Agertz et al.|(|2013f. Briefly, 
following Hemquist| ( |1993| ) and |Springel| ( |2000| l (see also |Springel| 
|et al.|2005| we cre^ a particle distribution representing a late type, 
star forming spiral galaxy embedded in an NFW dark matter halo 
( [Navarro et al.||T996t . The dark matter halo has a concentration 
parameter c = 10 and virial circular velocity, measured at the over¬ 
density 200pcrit, V 200 = 150kms“^, which translates to a halo 
virial mass M 200 = 1.1 x 10^^ Mq. The total baryonic disc mass 
is Mdisc = 4.5 X 10^° Mq with 20% in gas. The bulge-to-disc 
mass ratio is B/D = 0.1. We assume exponential profiles for 
the stellar and gaseous components and adopt a disc scale length 
Td = 3.4 kpc and scale height h — O.lrd for both. The bulge mass 
profile is that of |Hemquist| ( | 1 990^ with scale-length a = O.lrd. The 
halo and stellar disk are represented by 10® particles each, and the 
bulge consists of 10® particles. 

We initialize the gaseous disc analytically on the AMR grid 
assuming an exponential profile. The galaxy is embedded in a hot 
(T = 10® K), tenuous (n = 10“®cm“®) gas halo enriched to 


Z = 10-^Zq, while the disc has solar abundance. The minimum 
AMR cell size reached in the simulations is Ax = 9 pc. 

We carry out two simulations, one without any stellar feed¬ 
back, adopting a local star formation efficiency per free-fall time 
(see e.g. [Agertz et al.|2013^ ts = 1%, and one with stellar feed¬ 
back but with cff = 10%. In the former case, the effect of feedback 
is implicit in the choice of eg, and in the latter case this is achieved 
by efficient stellar feedback. Both simulations therefor have simi¬ 
lar star formation histories and normalizations of the Esfr — Egaa 
(Kennicutt-Schmidt) relation (see discussion injAg ertz et al.|20T^ 
[Agertz & Kjavtsov|2014), but reach this state in different ways. 


3 RESULTS 

3.1 Classical stability analysis 

The left hand side of Fig.[T] shows projected gas surface density 
maps of the two simulated galactic discs at t = 140 Myr. With¬ 
out stellar feedback, the disc violently fragments into a long-lived 
population of massive star forming GMCs, whereas the clouds are 
rapidly dispersed, and reformed, in the feedback regulated case. In 
the right hand panel we show radial profiles of the classical Toomre 
Q for gas and stars, as well as the two-component thin and thick 
disc stability parameter, Qthin and Qthick respectively (see § [2.1^ . 

We find all Q radial profiles to be almost indistinguishable be¬ 
tween the two models, and note that they are in agreement with de¬ 
rived values from well resolved spiral galaxies, see e.g. the detailed 
analysis b y [Leroy et al.[(|2008^ of th e THINGS galaxies (see also 
figure 5 in [Romeo & Wiegert|2011[ [Romeo & Falstad|2013) ; e.g. 
Qthick ~ 1-5 — 3 for r > 2 kpc, i.e. outside the bulge, indicative 
of gravitational stability to axisymmetric waves. The similar values 
of Q(r), all being significantly larger than unity, arise despite the 
morphological states of the discs being markedly different. 

[Leroy et al.j ( |2008[ ( concluded that the large observed val¬ 
ues of Q, indicative of local stabilit}]^ put doubts on the role of 
gravitational instabilities in alone controlling the local star forma¬ 
tion efficiency. However, as we have been alluding to in previous 
sections, our analysis is done on a fixed, and rather large, scale 
(io = 0.5 kpc) chosen to coincide with the spatial resolution of the 
THINGS survey. The fact that different indicators of stability fail 
to separate the simulated galaxies means that galactic dynamics, 
on the chosen smoothing scale, is the same. To probe differences 
further we need to adopt a scale dependent stability analysis. 


3.2 Accounting for turbulent scalings 

In this section we characterize the effect of turbulence on the grav¬ 
itational stability of the simulated discs directly via the dispersion 
relation for gas and stars (Eq.|^, as well as the relation for the cou¬ 
pled system (tu?. in Eq.[l2|. 

We estimate E(f) and a{£) for gas and stars from the sim¬ 
ulations at different simulation times as follows; assume the disc 
plane coincides with the a:i/-plane, i.e. the axis of rotation is paral¬ 
lel with the jz-axis. Centered in the mid-plane of the disc (z = 0), 
we place cubes in a lattice configuration in the range {x,y) G 
{ — 15,15} kpc, with a lattice spacing of Al — 100 pc. Starting 
from a cube size of £ = 36 pc, we incrementally increase this value 


® the cooling radius scales as rcooi ~ (Z/Zq + 0. 01)"*^'^^ pc 

for a supernova explosion with energy Egfq = 10®^ erg (e.g. [Cioffi et al.[ 
[1988[[Thomton et al.[199^[Kim & Ostriker[2014| 


^ Note that this refers to stability against axisymmetric perturbations and 
that stability against no«-axisymmetric perturbations is not guaranteed (e.g. 
[Griv & Gedalin[2012| 
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r [kpc] 

Figure 1. (Left) Projected gas density maps covering a region 36 X 36 kpc in size at f = 140 Myr of the simulation without (top) and with (bottom) stellar 
feedback. (Right) Mass weighted average profiles, computed in radial bins of size Ar = 0.5 kpc of the classical Toomre Q (= (Tk/ttGE) for gas and stars, 
the joint stability parameter Qthin and the joint parameter accounting for disc thickness Qthick (Romeo & Wiegert 2011; Romeo & Falstad 2013). 


by 2Ax = 36 pc up to £ = 0.5 kpc, after which we increase the 
cube sizes in steps of 50% of the current value up to f ~ 10 kpc, 
simply to reduce the computational cost. On each scale we compute 
Ei(f) and for each component (gas and stars) for every cube 
i. We then define S(^) and a{£) to be the averages of all cubes at 
any given scale, and consider these quantities to be representative 
of the typical surface density and velocity dispersion. We have con¬ 
firmed that sampling the discs in different way, e.g. randomly, has 
a negligible effect on the results presented below. 


From now on we focus our analysis on the average scale de¬ 
pendent characteristics centered in an azimuth defined by 4 kpc < 
r < 5 kpc. Furthermore, as the gaseous velocity dispersion is 
found not be isotropic (see also discussion in |Agertz et al.|2009a| ), 
we will from here on consider the radial velocity dispersion aaif) 
for both gas and stars, as this is the relevant component entering the 
dispersion relation for non-isotropic velocity ellipsoids. Note also 
that the sound speed of the gas Cs ^ an in a mass weighted sense; 
the cold ISM is dominated by super-sonic motions, and we can in 
principle omit Cs in the stability analysis. 


3.2.1 Density evolution and spatial scaling 

Fig. 0 shows the time evolution (t = 60,100 and 140 Myr) of 
Egas(f) and Estar(-^) for the fragmenting (no feedback) and feed¬ 
back regulated discs. Both simulations show “saturated” quantities 
above £ ~ 1 kpc, above which little variation with scale is found. 


confirming the strikingly similar Q{r) values found in §3.1 On 
smaller scales, the two models evolve differently; without feedback 
at f = 140 Myr (and as well for earlier times, but on scales of a few 
100 pc), the population of bound star forming clouds, arising from 
violent fragmentation, leads to a steeply decreasing Egas(£) with 
increasing £, while the feedback regulated case shows a rather flat 
Egaa(£). Quantified using the power-law relation inEq.[^(E oc £“), 
the clumpy galaxy approaches a ~ — 1 and the feedback regulated 
a ~ 0 for £ < 1 kpc. Significant scatter exists in the data, owing to 
the structured nature of the ISM and the fact that we do not bias the 
analysis to any specific regions of the disc0 or phases of the gas, 
other than the cold ISM. 

Eatar(£) is increasing at scales larger than ~ 100 pc. This is 


® For example, by not enforcing the analysis to be centered only on the 
dense cloud population in the fragmented disc, we do not measure a mono- 
tonically increasing S(£) as £ —> 0. 
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Figure 2. The scale dependence of the gas and stellar surface density S in the models adopting no feedback (top) and feedback (bottom) for, from left to right, 
t = 60,100 and 140 Myr. 
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Figure 3. The scale dependence of the gas and stellar radial velocity dispersion in the models adopting no feedback (top) and feedback (bottom) for, from 
left to right, t = 60,100 and 140 Myr. 
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Figure 4. The instability growth rate, accounting for scale dependent aiy.') and S(£), for gas (solid lines) and stars (dashed lines) for no feedback (top) and 
feedback (bottom). When turbulent scalings are accounted for in the analysis, the feedback simulation is found to be on average stable on all scales, with a 
Toomre-like behaviour (see main text), whereas the model lacking feedback is characterized as unstable on small scales (f < a few 100 pc), in accordance 
with the fragmented gas morphology with star forming clouds of sizes ~ 100 pc. 


due to the analysis being done in mid-plane centered cubes, for 
which scales £ ~ 1 kpc do not encapsulate all of the (kinematically 
hot) stars present in the disc. A similar conclusion may be drawn 
for the cold gas in the feedback regulated simulation, where the 
small scale turbulence, in an average sense, allows for a positive 
value of a. Adopting a scale dependent analysis for both stars and 
gas hence automatically accounts for disc thickness, without the 
need to adopt an effective surface density Sefi = S/(l + kh), as 
is done in Eq.[^ 


3.2.2 Velocity dispersion evolution and spatial scaling 

Fig. 0 shows the time evolution, t = 60,100 and 140 Myr, of 
o'fl,star(^) and crij,gas (£) for the fragmenting (no feedback) and 
feedback regulated discs. The stars show a roughly scale indepen¬ 
dent value of Or ~ 30kms“^ in both disc models. For the cold 
gas we find that the velocity dispersions are, after an initial tran¬ 
sient, quite similar both in magnitude and dependency on scale re¬ 
gardless of the presence of feedback or not. In both models, and 
for t > 60 Myr, (7_R,gaa(^) oc on scales £ < 300 — 400 pc 
in the feedback model, and I < 100 pc in the model without, in 
excellent agreement with the observed Larson-like scaling relevant 
for GMCs (see §[TJ. 

Still at times f > 60 Myr, and on large scales (a few 100 pc < 
I < 4 kpc), the slope of ( 7 _R,gas(-f) transitions from & ~ 1/2, into 
a flatter profile, with 6 ~ 1/5, in the case without feedback. In 


the feedback regulated galaxy we measure h ~ 1/2 up to almost 
~ 1 kpc, with a flattening on scales > 1 kpc. The slopes, and the 
scales over which the relation is measured in the feedback model, 
are in excellent agreement with observations of the Galactic HI 
linewidth-scale relation. For example, |Kim et JI] ( |2007| l studied 
this relation for individual HI clouds in the Large Magellanic Cloud 
(LMC) and found b~0.3 — 0.5on scales < 0.5 — 1 kpc (down to 
parsec scales). 

We emphasize that both models show values of cr ~ 10 — 
20 km on large scales, in agreement with both HI and CO ob¬ 
servations in local spiral galaxies ( [Tamburro et al.|[2009[ |Caldu-| 
|Primo et al.|2013l >, despite the different nature in ISM turbulence 
driving on small scales. In the case of feedback-driven turbulence, 
star forming clouds are rapidly destroyed by internal process and 
the gas is dispersed, whereas disordered gas motions are driven 
by gravitational instabilities and galactic shear, leading to clump- 
clump interactions (e.g. Jog^()strikeijl98^ |Gammie et al.|1991t 


[Agertz et al.|2009^ [Tasker & Tan|2009 1 in the case without feed¬ 

back. 


3.2.3 The dispersion relation 

Fig. 0 shows the resulting mass weighted average dispersion rela¬ 
tions, and the associated scatter, for the two galactic models for 
the stars cold gas (Wgas) and the coupled system (w?.). 

We remind the reader that the average aP(£) relations now self- 
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consistently account for the actual scale dependent values of a{(,) 
and E(^) present in each analyzed region of the galaxies. 

While featuring an almost identical Q{r) (for t = 140 Myr, 
& when analyzed on ~ 1 kpc scales, we can here identify 
the level of stability on all scales. Both discs are indeed stable or 
marginally stable > 0 ) on large scales, and have in fact almost 
identical relations for £> 1 kpc. 

The fragmented disc without stellar feedback always shows 
cUgas < 0 and uj1_ < 0 on small scales, as suspected simply via 
visual inspection of the top left panel of Fig.[T] At early times 
(t — 60 Myr), when fragmentation has just occured, the scale of 
instability is formally as large as ~ 2 kpc, with scales below a few 
~ 100 pc being the most unstable. At subsequent times, wi < 0 
on scales £ < 100 — 200 pc, roughly the maximum size of GMCs, 
or more correctly GMAs (Giant Molecular Associations) in this 
model. This is also the scale at which we measured a transition 
in the Larson-like scaling of the gas velocity dispersions, from 
b ~ 1/2 into b ~ 1/5 (§ [TT^ . 

The above conclusions are in stark contrast to the feedback 
regulated case which is, at least on average, stable or marginally 
stable on all scales. Note that individual patches of the disc on small 
scales can be unstable, leading to star formation, as is evident for 
^gas(^ ^ 100 pc) at a > 0.5 a level. 

This analysis underscores the necessity to extend a traditional 
stability analysis with scale dependent variables to account for the 
typical average velocity and density structures that exists in the 
ISM. Furthermore, many classical concepts, such as a well defined 
fastest growing mode for a Toomre unstable (Q < 1) disc, may no 
longer exist for the emerging scalings of E(^) and a{£) introduced 
by strong fragmentation, as pointed out by |Romeo et al.| ( [2010| l and 
[Romeo & AgertzH2014^ . 


3.3 Mapping out the stability regimes 

How do the results in the previous sections connect to observed 
turbulent scaling in the ISM and results from high resolution nu¬ 
merical simulations? Fig.|5] shows the stability map of turbulence 
( [Romeo et al.||20iQl l, where the axes denote the a and b power- 
law exponents for the surface density and the velocity dispersion 
(Eq.§. 

Following the analysis in [Romeo et al.| l |2010^ (see also [Hoff^ 
[mann & Romeo[2012[ l, we have identified, and indicated in the fig¬ 
ure, three regimes of stability: 

• Regime A: For b < | (1 -F a) and —2 < a < 1, the stability 
of the disc is controlled by Qo = ttcto /ttGEo (i.e. the Q measured 
on the fiducial smoothing scale £o, see § |2.2^ : the disc is stable at 
all scales if and only if Qo > Qq, where depends on a, b and 
£o, and must be derived from Eq.|^for a single component system. 
This is the domain of HI turbulence. Both HI observations and 
high-resolution simulations of supersonic turbulence are consistent 
with the scaling a ~ fo. In the case when a = b, the ratio a{£)/T,{£) 
is a constant and the local stability criterion degenerates into the 
classical Toomre case, Qo > 1, as if the disc was non-turbulent 
and infinitesimally thin. 

• Regime C: For 6 > | (1 -F a) and — 2 < o < 1, the stability of 
the disc is no longer controlled by Qo'- the disc is always unstable 
on small scales (i.e. as / —>■ 0) and stable on large scales. 

• Regime B: For b = | (1 + a) and —2 < a < 1, the disc 
is in a phase of transition between Toomre-like stability (Regime 
A) and instability om small scales (Regime C). This is the domain 
typically observed for molecular gas in GMCs. 



a 


Figure 5. The stability map of turbulence, with typical values derived from 
observations and numerical simulations, see the main text for a comprehen¬ 
sive discussion. The blue, red and black points indicate the regimes found 
in the simulated disc galaxies on both small and large scales. On large 
scales (f > 1 kpc) the feedback regulated disc and the fragmented disc 
converge to roughly the same average S and a, and both well described by 
a Toomre-like stability criterion (Regime A, see text). Turbulent scaling on 
small scales pushes the fragmented disc into a regime typical observed for 
GMCs (Regime B), and a regime where the classical Toomre analysis is no 
longer valid (Regime C). The disc with feedback-driven turbulence is, in 
a statistical sense, stabilized on small scales, and shows scalings in broad 
agreement with HI and CO observations. 


A few well studied examples from theory/simulations and ob¬ 
servations in the literature are shown in the figure: 


• (L), (a, 6 ) = (—0.1,0.38): The original scaling relations in 
GMCs found by [Larson'| ( [1981[ l. 

• (S), (a, b) = (0,0.50 zF 0.05): The scaling relations in GMCs 
found by [Solomon et al.|jl987[ l. These values are what is usually 
referred to as “Larson’s scaling laws”. 

• (K), (a, 6) = (|i I)' The result of high-resolution simula¬ 
tions of supersonic turbulence ( [Kritsuk et al.|201^ Kritsuk, private 
communication). 

• (F), (a, b) = (0.44 ± 0.14, 0.49 ± 0.02): A prediction based 
on state-of-the-art simulations of super-sonic turbulence with com¬ 
pressive driving ( [Federrath||2013| Federrath, private communica¬ 
tion). For solenoidal driving, (a, b) = (0.58 zF 0.03,0.48 zF 0.02). 


• (E), (a, b) = (—1, 0.5): Investigated by Elmegreen 


• (HI): Typical range of values derived from observed HI inten¬ 
sity fluctuations in disc galaxies ([Lazarian & Pogosyan|2000|[Kim[ 
[et al.|2007[[Roy et al.|2008) . 


A number of typical (a, b) pairs measured in the two simu¬ 
lations are indicated in Eig.[^ In the case of no feedback, lead¬ 
ing to strong fragmentation, the disc features a ~ —1 to 0 , and 
6 ~ 1/5 — 1/2 on small/intermediate scales. This puts the disc in 
regime B and C, with the latter meaning that the disc is always un¬ 
stable on small scales. In fact, by studying the dispersion relations 
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(Wgaa(^) ora;?.(^))inFig. 0 no well defined minimum exists, with 
the smallest numerically resolved scale, here on the order of a cell 
size ~ Ax, being the fastest growing mode. On large scales, we 
measure a transition into Regime A with (a, b) ~ (0, 0), i.e. the 
classical Toomre case with a well defined average surface density 
and velocity dispersion. 

The feedback driven simulation shows the same behaviour on 
large scales > 1 kpc) as the fragmented counterpart, coincid¬ 
ing with the classical Toomre case. On small scales < 1 kpc), 
the simulations diverge, as discussed already in § |3.2.3| Here feed¬ 
back creates an, on average, marginally stable ISM, still in regime 
A, with a and b values compatible with observations of the cold 
galactic ISM (HI and CO). Using these power law exponents, one 
can define a Q stability parameter from Eq.|^or[T^that quantifies 
the threshold for instability, in contrast to the case where the disc 
fragments due to the lack of stellar feedback support. 


4 DISCUSSION AND CONCLUSIONS 


Gravitational instabilities are thought to play an important role in 
galaxy evolution and star formation. In this work we have investi¬ 
gated, using numerical simulations of Milky Way-like disc galax¬ 
ies, the nature of turbulence in the ISM and how this affects the 
gravitational stability of galaxies. Clumpy/turbulent discs are dy¬ 
namically similar to gas discs with scale-dependent surface densi¬ 
ties and velocity dispersions, i.e. E(£) ~ l°‘ and ~ re¬ 
spectively, where i is the physical scale. By taking these “turbu¬ 
lent” scaling relations into account in the disc stability analysis, a 
wide variety of non-classical stability properties arise (|Romeo et al.| 

|2m^ . 

In order to quantify the kind of turbulent scaling that de¬ 
velop in the ISM of disc galaxies we used numerical simulations 
of multi-component. Milky Way-like, galactic discs. We studied 
two markedly different evolutionary scenarios using the same ini¬ 
tial condition: ( 1 ) no stellar feedback, leading to complete disc frag¬ 
mentation and (2) efficient stellar feedback leading to driven ISM 
turbulence. By accounting for the measured turbulent scalings in 
the stability analysis, we could more robustly characterize the level 
stability as a function of physical scale. Our key results are sum¬ 
marized below: 


(i) Our two different models of galaxy evolution lead to discs 
with nearly identical stability properties when quantified on kpc- 
scales by the classical Toomre Q for stars and gas separately 
{Q S> 1), as well as jointly (Qstars+gas > 1, accounting for disc 
thickness); the discs are stable at all radii. This notion is in good 
agreement with observations of nearby disc galaxies, e.g. THINGS 
( [Leroy et al.|2008l l where all discs, smoothed on kpc scales, fea¬ 
tured values of Qgaa and Q* well above unity (see also [Romeo &| 
|Wiegert|2011| (, raising doubts on the role of gravitational instabili¬ 
ties in star formation. 

(ii) The two models feature markedly different average E(£) on 

scales £ < Ikpc, with a steepening of E(^) into ~ £~^ in the 
fragmented disc and ~ £0-i/3 feedback regulated disc. 

Less of a diffe rence was found for a'(£), with a La rson-like scaling 


Larson 


1981 


Solomon et al. 


1987 I in both models. 


(a{£) ~ 

(hi) Introducing scale dependent variables in the multi- 
component stability analysis leads to a more robust characteriza¬ 
tion of the level of instability. The feedback driven model is, on 
average, stable or marginally stable on all scales, in contrast to 
the model without feedback, for which we can clearly identify the 
scales (£ < 100 pc) where gravitational instability, leading to cloud 


formation, typically occurs. Large scales (^ > 1 kpc) show almost 
identical stability properties in both models, as the surface den¬ 
sity and velocity dispersion “saturate” into more well-defined, non- 
turbulent quantities, explaining the similarity between the markedly 
different ISM models when adopting a traditional large scale stabil¬ 
ity analysis. 

(iv) The disc stabilized by stellar feedback can still be, in a sta¬ 
tistical sense, well described by a Toomre-like Q stability thresh¬ 
old. This is no longer true for the violently fragmenting disc, which 
enters a regime where a Q-like parameter loses its meaning and 
small scales are asymptotically more and more unstable, and sta¬ 
bility can only occur on large scales (see alsojRomeo et al.|2010). 


[Hopkins & Christiansenj ( [2013^ investigated gravitational 
fragmentation in turbulent media, arguing that turbulent flows are 
always unstable on some scale, given enough time, as a broad spec¬ 
trum of stochastic density fluctuations exists that can produce rare, 
but gravitationally unstable regions. This notion is compatible with 
our results; by defining the typical densities and velocities existing 
at some scale in a mass weighted sense, we are biasing ourselves 
towards dense regions. This gives us the typical structures, by mass, 
that exists in the ISM. However, the scatter in the dispersion rela¬ 
tion for the model with feedback-driven turbulence (Fig.j^ can be 
substantial on small scales at any time, due to the wide range of 
densities in supersonically turbulent flows (e.g. [Padoan et al.[1997[ 
[Mac Low & Klessen[2004[[Kritsuk et al.[2007^ . This means parts of 
the disc will be unstable, although the corresponding mass in that 
component may not necessarily be dominant. 

In this work we have only considered Milky Way-like galax¬ 
ies, with gas fractions typical of ^ = 0 discs. In high redshift 
counterparts, with gas fractions several times greater jTacconi et al.[ 
|20T0|, turbulence may play an even greater role in shaping the 
ISM. Indeed, the observed high levels of turbulence (e.g. [Forster [ 
[Schreiber et al.[[2009[( is thought to be responsible for the highly 


clumpy morphologies observed (e.g. [Elmegreen et al.[[2009 

[Tac- 

coni et al. [2010 Agertzetal. 

2009b[[Dekel et al.[2009| Genze 

et al. 

201 1| Romeo & Agertz[20T 

I. Current astronomical facilities such 


as ALMA can resolve the scaling properties of galactic turbulence 
in the cold molecular gas of high redshift systems, hence reveal¬ 
ing the interplay between gravitational instability and turbulence in 
more extreme environments. 

It is important to consider the limitations and assumptions be¬ 
hind the analysis presented in this work. The dispersion relation 
(Eq.j^andjT} is formally only describing stability against local ax- 
isymmetric perturbations, i.e. it assumes that kR 3> 1. This is the 
short-wavelength approximation and is satisfied in our work as we 
only carry out our analysis at i? = 4 kpc, and with k = 2tt/£ this 
holds on all scales we have considered. However, the assumption 
of axisymmetry is not true in general, and nonaxisymmetric pertur¬ 
bations are thought to have a greater destabilizing effect, i.e. discs 
that are formally stable (Q > 1) can be locally unstable against 
such perturbations. Local nonaxisymmetric stability criteria are far 
more complex than Toomre’s criterion as they depend critically on 
how tightly wound the perturbations are, and any such criteria can¬ 
not, in general, be expressed in terms of a single parameter akin to 
Toomre’s Q (e.g.[Lau & Bertin[1978[[Bertin et al.[1989[[jog[1992[ 
[Griv & Gedalin|2012^ . We leave an analysis accounting for more 

general perturbations for future work. 

In future work (Grisdale et al. in prep) we will extend the anal¬ 
ysis presented in this paper to quantify the gravitational stability for 
different gas phases, regions of the disc, as a function of local tur- 
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bulence driving strength etc., and how this connects to properties 
of observed galaxies. 
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